# -*- coding: utf-8 -*-
# Code sous licence CeCill
# http://www.cecill.info/licences.en.html
# This is a license compatible with the GNU GPL license
# and that conforms to the principles of French law.
# Authors:
#       Nicolas Bruyant, CNRS,LNCMI
#       Benjamin Edwards, CNRS,LNCMI
"""
Import Filter
=============

Utility functions to import datas and calculate field from pick-up datas

Functions index:
    loadfile
    field_integration
    field_adjust

Todo:
----
    Improve loadfile function to allow column naming

"""
__version__ = "Version 0.1, SVN-$Revision: 18 $"
"""
Revision allow a good tracking of data processing  errors
Refers to SVN revisions: DO NOT edit manually
"""

import numpy as np

from scipy import integrate,mean,arange,polyfit



def loadfileFFC(filename,sample_rate=500E3,nb_col=5):
    """
    Load raw data file from a file of type binary FFC

    Parameters
    ==========
        Filename,
        Sample rate of the datas,
        Number of column of the datafile

    Returns
    =======
        Numpy array containing datas

    Examples
    ========
    >>> inputfile='/Users/nico/Documents/Manips/Superadiance/Documents/Toulouse experiment/Data/jva004'
    >>> temp=loadfile(inputfile,10E6,4)

    """
    with open(str(filename),'rb') as fd:
        temp=np.fromfile(fd,np.float32,-1,"")
    temp.resize(len(temp)/nb_col,nb_col)
    return temp


def loadfileGage(filename, sample_rate=1E6):
    """
    Load raw data file from a file of type binary GageScope

    Parameters
    ==========
        Filename,
        Sample rate of the datas,
        Number of column of the datafile

    Returns
    =======
        Numpy array containing datas

    Examples
    ========
    >>> filename="/Users/Nicolas/Documents/Projets/TDO/post traitement/tracking freq/s120725d.sig"
    >>> temp=loadfileGage(filename)

    """
    with open(str(filename), 'rb') as fd:
        fd.seek(512)   # Skip useless header
        #TODO: extract useful datas ie samplerate
        temp = np.fromfile(fd, np.int16, -1, "")
    return temp


def field_integration(dbdt_data,sample_rate=500E3,method='all',
                      offset_calculation_length=100):
    """
    Integrate pick-up signal for the field

    Parameters
    ==========
        dbdt datas,
        Sample rate,
        Method of offset substraction (Optional):
            -"all": Mean of the whole dataset (default)
            -"start": Mean of the first N samples
            -"end": Mean of the last N samples
        offset_calculation_length: number of sample for offset calculation

    Returns
    ======
        Integrated pick-up signal

    """
    if method=='end':
        offset=np.mean(dbdt_data[-offset_calculation_length:])
    elif method=='start':
        offset=np.mean(dbdt_data[0:offset_calculation_length])
    else:
        offset=np.mean(dbdt_data)

    temps=np.arange(0,len(dbdt_data))*1/sample_rate
    field=integrate.cumtrapz(dbdt_data-offset,temps)
    return field

def field_adjust(champ):
    """
    Adjust the calculated field signal to remove residual offset

    Parameters
    ==========
        champ:  array
            Integrated pick-up signal

    Returns
    =======
        Corrected signal

    TODO
    ====
    utiliser le fit automatique sur toute la plage
    """
    # threshold pour un champ superieur à 10% du max
    threshold=max(champ)/20

    #on se place 10% ech avant le threshold
    for i in range(len(champ)):
       if champ[i]>threshold:
            break
    else:
       print "pb de seuil"
    imin=int(i*0.9)
    for i in range(i,len(champ)):
        if champ[i]<threshold:
            break
    else:
        print "pb de seuil"

    imax=int(i*1.1)
    m1=np.zeros(2)
    m2=np.zeros(2)
    # on fitte en lineaire
    m1[1]=mean(champ[:imin]/arange(0,imin))
    m2[1]=mean(champ[imax:]/arange(imax, len(champ)))
    m1[0]=0
    m2[0]=0

#    Marche pas ...
#    x=concatenate(arange(0,imin),arange(imax, len(champ)))
#    y=champ[:imin],champ[imax:]
#    m=polyfit(x,y,1)
#


    #on choisis la "meilleure", i.e. la plus longue plage de fit
    if imin > (len(champ)-imax):
        champ_corr=champ-(m1[0]+m1[1]*arange(0,len(champ)))
        print m1
    else:
        champ_corr=champ-(m2[0]+m2[1]*arange(0,len(champ)))
        print m2
    print "max field= %e T" % max(champ_corr)
    return champ_corr
